/*************************************************************************
Cephes Math Library Release 2.8:  June, 2000
Copyright by Stephen L. Moshier

Contributors:
    * Sergey Bochkanov (ALGLIB project). Translation from C to
      pseudocode.

See subroutines comments for additional copyrights.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

- Redistributions of source code must retain the above copyright
  notice, this list of conditions and the following disclaimer.

- Redistributions in binary form must reproduce the above copyright
  notice, this list of conditions and the following disclaimer listed
  in this license in the documentation and/or other materials
  provided with the distribution.

- Neither the name of the copyright holders nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*************************************************************************/

#include <stdafx.h>
#include "bessel.h"

static void besselmfirstcheb(double c, double& b0, double& b1, double& b2);
static void besselmnextcheb(double x,
     double c,
     double& b0,
     double& b1,
     double& b2);
static void besselm1firstcheb(double c, double& b0, double& b1, double& b2);
static void besselm1nextcheb(double x,
     double c,
     double& b0,
     double& b1,
     double& b2);
static void besselasympt0(double x, double& pzero, double& qzero);
static void besselasympt1(double x, double& pzero, double& qzero);

/*************************************************************************
Bessel function of order zero

Returns Bessel function of order zero of the argument.

The domain is divided into the intervals [0, 5] and
(5, infinity). In the first interval the following rational
approximation is used:


       2         2
(w - r  ) (w - r  ) P (w) / Q (w)
      1         2    3       8

           2
where w = x  and the two r's are zeros of the function.

In the second interval, the Hankel asymptotic expansion
is employed with two rational functions of degree 6/6
and 7/7.

ACCURACY:

                     Absolute error:
arithmetic   domain     # trials      peak         rms
   IEEE      0, 30       60000       4.2e-16     1.1e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*************************************************************************/
double besselj0(double x)
{
    double result;
    double xsq;
    double nn;
    double pzero;
    double qzero;
    double p1;
    double q1;

    if( x<0 )
    {
        x = -x;
    }
    if( x>8.0 )
    {
        besselasympt0(x, pzero, qzero);
        nn = x-ap::pi()/4;
        result = sqrt(2/ap::pi()/x)*(pzero*cos(nn)-qzero*sin(nn));
        return result;
    }
    xsq = ap::sqr(x);
    p1 = 26857.86856980014981415848441;
    p1 = -40504123.71833132706360663322+xsq*p1;
    p1 = 25071582855.36881945555156435+xsq*p1;
    p1 = -8085222034853.793871199468171+xsq*p1;
    p1 = 1434354939140344.111664316553+xsq*p1;
    p1 = -136762035308817138.6865416609+xsq*p1;
    p1 = 6382059341072356562.289432465+xsq*p1;
    p1 = -117915762910761053603.8440800+xsq*p1;
    p1 = 493378725179413356181.6813446+xsq*p1;
    q1 = 1.0;
    q1 = 1363.063652328970604442810507+xsq*q1;
    q1 = 1114636.098462985378182402543+xsq*q1;
    q1 = 669998767.2982239671814028660+xsq*q1;
    q1 = 312304311494.1213172572469442+xsq*q1;
    q1 = 112775673967979.8507056031594+xsq*q1;
    q1 = 30246356167094626.98627330784+xsq*q1;
    q1 = 5428918384092285160.200195092+xsq*q1;
    q1 = 493378725179413356211.3278438+xsq*q1;
    result = p1/q1;
    return result;
}


/*************************************************************************
Bessel function of order one

Returns Bessel function of order one of the argument.

The domain is divided into the intervals [0, 8] and
(8, infinity). In the first interval a 24 term Chebyshev
expansion is used. In the second, the asymptotic
trigonometric representation is employed using two
rational functions of degree 5/5.

ACCURACY:

                     Absolute error:
arithmetic   domain      # trials      peak         rms
   IEEE      0, 30       30000       2.6e-16     1.1e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*************************************************************************/
double besselj1(double x)
{
    double result;
    double s;
    double xsq;
    double nn;
    double pzero;
    double qzero;
    double p1;
    double q1;

    s = ap::sign(x);
    if( x<0 )
    {
        x = -x;
    }
    if( x>8.0 )
    {
        besselasympt1(x, pzero, qzero);
        nn = x-3*ap::pi()/4;
        result = sqrt(2/ap::pi()/x)*(pzero*cos(nn)-qzero*sin(nn));
        if( s<0 )
        {
            result = -result;
        }
        return result;
    }
    xsq = ap::sqr(x);
    p1 = 2701.122710892323414856790990;
    p1 = -4695753.530642995859767162166+xsq*p1;
    p1 = 3413234182.301700539091292655+xsq*p1;
    p1 = -1322983480332.126453125473247+xsq*p1;
    p1 = 290879526383477.5409737601689+xsq*p1;
    p1 = -35888175699101060.50743641413+xsq*p1;
    p1 = 2316433580634002297.931815435+xsq*p1;
    p1 = -66721065689249162980.20941484+xsq*p1;
    p1 = 581199354001606143928.050809+xsq*p1;
    q1 = 1.0;
    q1 = 1606.931573481487801970916749+xsq*q1;
    q1 = 1501793.594998585505921097578+xsq*q1;
    q1 = 1013863514.358673989967045588+xsq*q1;
    q1 = 524371026216.7649715406728642+xsq*q1;
    q1 = 208166122130760.7351240184229+xsq*q1;
    q1 = 60920613989175217.46105196863+xsq*q1;
    q1 = 11857707121903209998.37113348+xsq*q1;
    q1 = 1162398708003212287858.529400+xsq*q1;
    result = s*x*p1/q1;
    return result;
}


/*************************************************************************
Bessel function of integer order

Returns Bessel function of order n, where n is a
(possibly negative) integer.

The ratio of jn(x) to j0(x) is computed by backward
recurrence.  First the ratio jn/jn-1 is found by a
continued fraction expansion.  Then the recurrence
relating successive orders is applied until j0 or j1 is
reached.

If n = 0 or 1 the routine for j0 or j1 is called
directly.

ACCURACY:

                     Absolute error:
arithmetic   range      # trials      peak         rms
   IEEE      0, 30        5000       4.4e-16     7.9e-17


Not suitable for large n or x. Use jv() (fractional order) instead.

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
double besseljn(int n, double x)
{
    double result;
    double pkm2;
    double pkm1;
    double pk;
    double xk;
    double r;
    double ans;
    int k;
    int sg;

    if( n<0 )
    {
        n = -n;
        if( n%2==0 )
        {
            sg = 1;
        }
        else
        {
            sg = -1;
        }
    }
    else
    {
        sg = 1;
    }
    if( x<0 )
    {
        if( n%2!=0 )
        {
            sg = -sg;
        }
        x = -x;
    }
    if( n==0 )
    {
        result = sg*besselj0(x);
        return result;
    }
    if( n==1 )
    {
        result = sg*besselj1(x);
        return result;
    }
    if( n==2 )
    {
        if( x==0 )
        {
            result = 0;
        }
        else
        {
            result = sg*(2.0*besselj1(x)/x-besselj0(x));
        }
        return result;
    }
    if( x<ap::machineepsilon )
    {
        result = 0;
        return result;
    }
    k = 53;
    pk = 2*(n+k);
    ans = pk;
    xk = x*x;
    do
    {
        pk = pk-2.0;
        ans = pk-xk/ans;
        k = k-1;
    }
    while(k!=0);
    ans = x/ans;
    pk = 1.0;
    pkm1 = 1.0/ans;
    k = n-1;
    r = 2*k;
    do
    {
        pkm2 = (pkm1*r-pk*x)/x;
        pk = pkm1;
        pkm1 = pkm2;
        r = r-2.0;
        k = k-1;
    }
    while(k!=0);
    if( fabs(pk)>fabs(pkm1) )
    {
        ans = besselj1(x)/pk;
    }
    else
    {
        ans = besselj0(x)/pkm1;
    }
    result = sg*ans;
    return result;
}


/*************************************************************************
Bessel function of the second kind, order zero

Returns Bessel function of the second kind, of order
zero, of the argument.

The domain is divided into the intervals [0, 5] and
(5, infinity). In the first interval a rational approximation
R(x) is employed to compute
  y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
Thus a call to j0() is required.

In the second interval, the Hankel asymptotic expansion
is employed with two rational functions of degree 6/6
and 7/7.



ACCURACY:

 Absolute error, when y0(x) < 1; else relative error:

arithmetic   domain     # trials      peak         rms
   IEEE      0, 30       30000       1.3e-15     1.6e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*************************************************************************/
double bessely0(double x)
{
    double result;
    double nn;
    double xsq;
    double pzero;
    double qzero;
    double p4;
    double q4;

    if( x>8.0 )
    {
        besselasympt0(x, pzero, qzero);
        nn = x-ap::pi()/4;
        result = sqrt(2/ap::pi()/x)*(pzero*sin(nn)+qzero*cos(nn));
        return result;
    }
    xsq = ap::sqr(x);
    p4 = -41370.35497933148554125235152;
    p4 = 59152134.65686889654273830069+xsq*p4;
    p4 = -34363712229.79040378171030138+xsq*p4;
    p4 = 10255208596863.94284509167421+xsq*p4;
    p4 = -1648605817185729.473122082537+xsq*p4;
    p4 = 137562431639934407.8571335453+xsq*p4;
    p4 = -5247065581112764941.297350814+xsq*p4;
    p4 = 65874732757195549259.99402049+xsq*p4;
    p4 = -27502866786291095837.01933175+xsq*p4;
    q4 = 1.0;
    q4 = 1282.452772478993804176329391+xsq*q4;
    q4 = 1001702.641288906265666651753+xsq*q4;
    q4 = 579512264.0700729537480087915+xsq*q4;
    q4 = 261306575504.1081249568482092+xsq*q4;
    q4 = 91620380340751.85262489147968+xsq*q4;
    q4 = 23928830434997818.57439356652+xsq*q4;
    q4 = 4192417043410839973.904769661+xsq*q4;
    q4 = 372645883898616588198.9980+xsq*q4;
    result = p4/q4+2/ap::pi()*besselj0(x)*log(x);
    return result;
}


/*************************************************************************
Bessel function of second kind of order one

Returns Bessel function of the second kind of order one
of the argument.

The domain is divided into the intervals [0, 8] and
(8, infinity). In the first interval a 25 term Chebyshev
expansion is used, and a call to j1() is required.
In the second, the asymptotic trigonometric representation
is employed using two rational functions of degree 5/5.

ACCURACY:

                     Absolute error:
arithmetic   domain      # trials      peak         rms
   IEEE      0, 30       30000       1.0e-15     1.3e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
*************************************************************************/
double bessely1(double x)
{
    double result;
    double nn;
    double xsq;
    double pzero;
    double qzero;
    double p4;
    double q4;

    if( x>8.0 )
    {
        besselasympt1(x, pzero, qzero);
        nn = x-3*ap::pi()/4;
        result = sqrt(2/ap::pi()/x)*(pzero*sin(nn)+qzero*cos(nn));
        return result;
    }
    xsq = ap::sqr(x);
    p4 = -2108847.540133123652824139923;
    p4 = 3639488548.124002058278999428+xsq*p4;
    p4 = -2580681702194.450950541426399+xsq*p4;
    p4 = 956993023992168.3481121552788+xsq*p4;
    p4 = -196588746272214065.8820322248+xsq*p4;
    p4 = 21931073399177975921.11427556+xsq*p4;
    p4 = -1212297555414509577913.561535+xsq*p4;
    p4 = 26554738314348543268942.48968+xsq*p4;
    p4 = -99637534243069222259967.44354+xsq*p4;
    q4 = 1.0;
    q4 = 1612.361029677000859332072312+xsq*q4;
    q4 = 1563282.754899580604737366452+xsq*q4;
    q4 = 1128686837.169442121732366891+xsq*q4;
    q4 = 646534088126.5275571961681500+xsq*q4;
    q4 = 297663212564727.6729292742282+xsq*q4;
    q4 = 108225825940881955.2553850180+xsq*q4;
    q4 = 29549879358971486742.90758119+xsq*q4;
    q4 = 5435310377188854170800.653097+xsq*q4;
    q4 = 508206736694124324531442.4152+xsq*q4;
    result = x*p4/q4+2/ap::pi()*(besselj1(x)*log(x)-1/x);
    return result;
}


/*************************************************************************
Bessel function of second kind of integer order

Returns Bessel function of order n, where n is a
(possibly negative) integer.

The function is evaluated by forward recurrence on
n, starting with values computed by the routines
y0() and y1().

If n = 0 or 1 the routine for y0 or y1 is called
directly.

ACCURACY:
                     Absolute error, except relative
                     when y > 1:
arithmetic   domain     # trials      peak         rms
   IEEE      0, 30       30000       3.4e-15     4.3e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
double besselyn(int n, double x)
{
    double result;
    int i;
    double a;
    double b;
    double tmp;
    double s;

    s = 1;
    if( n<0 )
    {
        n = -n;
        if( n%2!=0 )
        {
            s = -1;
        }
    }
    if( n==0 )
    {
        result = bessely0(x);
        return result;
    }
    if( n==1 )
    {
        result = s*bessely1(x);
        return result;
    }
    a = bessely0(x);
    b = bessely1(x);
    for(i = 1; i <= n-1; i++)
    {
        tmp = b;
        b = 2*i/x*b-a;
        a = tmp;
    }
    result = s*b;
    return result;
}


/*************************************************************************
Modified Bessel function of order zero

Returns modified Bessel function of order zero of the
argument.

The function is defined as i0(x) = j0( ix ).

The range is partitioned into the two intervals [0,8] and
(8, infinity).  Chebyshev polynomial expansions are employed
in each interval.

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      0,30        30000       5.8e-16     1.4e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
double besseli0(double x)
{
    double result;
    double y;
    double v;
    double z;
    double b0;
    double b1;
    double b2;

    if( x<0 )
    {
        x = -x;
    }
    if( x<=8.0 )
    {
        y = x/2.0-2.0;
        besselmfirstcheb(-4.41534164647933937950E-18, b0, b1, b2);
        besselmnextcheb(y, 3.33079451882223809783E-17, b0, b1, b2);
        besselmnextcheb(y, -2.43127984654795469359E-16, b0, b1, b2);
        besselmnextcheb(y, 1.71539128555513303061E-15, b0, b1, b2);
        besselmnextcheb(y, -1.16853328779934516808E-14, b0, b1, b2);
        besselmnextcheb(y, 7.67618549860493561688E-14, b0, b1, b2);
        besselmnextcheb(y, -4.85644678311192946090E-13, b0, b1, b2);
        besselmnextcheb(y, 2.95505266312963983461E-12, b0, b1, b2);
        besselmnextcheb(y, -1.72682629144155570723E-11, b0, b1, b2);
        besselmnextcheb(y, 9.67580903537323691224E-11, b0, b1, b2);
        besselmnextcheb(y, -5.18979560163526290666E-10, b0, b1, b2);
        besselmnextcheb(y, 2.65982372468238665035E-9, b0, b1, b2);
        besselmnextcheb(y, -1.30002500998624804212E-8, b0, b1, b2);
        besselmnextcheb(y, 6.04699502254191894932E-8, b0, b1, b2);
        besselmnextcheb(y, -2.67079385394061173391E-7, b0, b1, b2);
        besselmnextcheb(y, 1.11738753912010371815E-6, b0, b1, b2);
        besselmnextcheb(y, -4.41673835845875056359E-6, b0, b1, b2);
        besselmnextcheb(y, 1.64484480707288970893E-5, b0, b1, b2);
        besselmnextcheb(y, -5.75419501008210370398E-5, b0, b1, b2);
        besselmnextcheb(y, 1.88502885095841655729E-4, b0, b1, b2);
        besselmnextcheb(y, -5.76375574538582365885E-4, b0, b1, b2);
        besselmnextcheb(y, 1.63947561694133579842E-3, b0, b1, b2);
        besselmnextcheb(y, -4.32430999505057594430E-3, b0, b1, b2);
        besselmnextcheb(y, 1.05464603945949983183E-2, b0, b1, b2);
        besselmnextcheb(y, -2.37374148058994688156E-2, b0, b1, b2);
        besselmnextcheb(y, 4.93052842396707084878E-2, b0, b1, b2);
        besselmnextcheb(y, -9.49010970480476444210E-2, b0, b1, b2);
        besselmnextcheb(y, 1.71620901522208775349E-1, b0, b1, b2);
        besselmnextcheb(y, -3.04682672343198398683E-1, b0, b1, b2);
        besselmnextcheb(y, 6.76795274409476084995E-1, b0, b1, b2);
        v = 0.5*(b0-b2);
        result = exp(x)*v;
        return result;
    }
    z = 32.0/x-2.0;
    besselmfirstcheb(-7.23318048787475395456E-18, b0, b1, b2);
    besselmnextcheb(z, -4.83050448594418207126E-18, b0, b1, b2);
    besselmnextcheb(z, 4.46562142029675999901E-17, b0, b1, b2);
    besselmnextcheb(z, 3.46122286769746109310E-17, b0, b1, b2);
    besselmnextcheb(z, -2.82762398051658348494E-16, b0, b1, b2);
    besselmnextcheb(z, -3.42548561967721913462E-16, b0, b1, b2);
    besselmnextcheb(z, 1.77256013305652638360E-15, b0, b1, b2);
    besselmnextcheb(z, 3.81168066935262242075E-15, b0, b1, b2);
    besselmnextcheb(z, -9.55484669882830764870E-15, b0, b1, b2);
    besselmnextcheb(z, -4.15056934728722208663E-14, b0, b1, b2);
    besselmnextcheb(z, 1.54008621752140982691E-14, b0, b1, b2);
    besselmnextcheb(z, 3.85277838274214270114E-13, b0, b1, b2);
    besselmnextcheb(z, 7.18012445138366623367E-13, b0, b1, b2);
    besselmnextcheb(z, -1.79417853150680611778E-12, b0, b1, b2);
    besselmnextcheb(z, -1.32158118404477131188E-11, b0, b1, b2);
    besselmnextcheb(z, -3.14991652796324136454E-11, b0, b1, b2);
    besselmnextcheb(z, 1.18891471078464383424E-11, b0, b1, b2);
    besselmnextcheb(z, 4.94060238822496958910E-10, b0, b1, b2);
    besselmnextcheb(z, 3.39623202570838634515E-9, b0, b1, b2);
    besselmnextcheb(z, 2.26666899049817806459E-8, b0, b1, b2);
    besselmnextcheb(z, 2.04891858946906374183E-7, b0, b1, b2);
    besselmnextcheb(z, 2.89137052083475648297E-6, b0, b1, b2);
    besselmnextcheb(z, 6.88975834691682398426E-5, b0, b1, b2);
    besselmnextcheb(z, 3.36911647825569408990E-3, b0, b1, b2);
    besselmnextcheb(z, 8.04490411014108831608E-1, b0, b1, b2);
    v = 0.5*(b0-b2);
    result = exp(x)*v/sqrt(x);
    return result;
}


/*************************************************************************
Modified Bessel function of order one

Returns modified Bessel function of order one of the
argument.

The function is defined as i1(x) = -i j1( ix ).

The range is partitioned into the two intervals [0,8] and
(8, infinity).  Chebyshev polynomial expansions are employed
in each interval.

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      0, 30       30000       1.9e-15     2.1e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1985, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
double besseli1(double x)
{
    double result;
    double y;
    double z;
    double v;
    double b0;
    double b1;
    double b2;

    z = fabs(x);
    if( z<=8.0 )
    {
        y = z/2.0-2.0;
        besselm1firstcheb(2.77791411276104639959E-18, b0, b1, b2);
        besselm1nextcheb(y, -2.11142121435816608115E-17, b0, b1, b2);
        besselm1nextcheb(y, 1.55363195773620046921E-16, b0, b1, b2);
        besselm1nextcheb(y, -1.10559694773538630805E-15, b0, b1, b2);
        besselm1nextcheb(y, 7.60068429473540693410E-15, b0, b1, b2);
        besselm1nextcheb(y, -5.04218550472791168711E-14, b0, b1, b2);
        besselm1nextcheb(y, 3.22379336594557470981E-13, b0, b1, b2);
        besselm1nextcheb(y, -1.98397439776494371520E-12, b0, b1, b2);
        besselm1nextcheb(y, 1.17361862988909016308E-11, b0, b1, b2);
        besselm1nextcheb(y, -6.66348972350202774223E-11, b0, b1, b2);
        besselm1nextcheb(y, 3.62559028155211703701E-10, b0, b1, b2);
        besselm1nextcheb(y, -1.88724975172282928790E-9, b0, b1, b2);
        besselm1nextcheb(y, 9.38153738649577178388E-9, b0, b1, b2);
        besselm1nextcheb(y, -4.44505912879632808065E-8, b0, b1, b2);
        besselm1nextcheb(y, 2.00329475355213526229E-7, b0, b1, b2);
        besselm1nextcheb(y, -8.56872026469545474066E-7, b0, b1, b2);
        besselm1nextcheb(y, 3.47025130813767847674E-6, b0, b1, b2);
        besselm1nextcheb(y, -1.32731636560394358279E-5, b0, b1, b2);
        besselm1nextcheb(y, 4.78156510755005422638E-5, b0, b1, b2);
        besselm1nextcheb(y, -1.61760815825896745588E-4, b0, b1, b2);
        besselm1nextcheb(y, 5.12285956168575772895E-4, b0, b1, b2);
        besselm1nextcheb(y, -1.51357245063125314899E-3, b0, b1, b2);
        besselm1nextcheb(y, 4.15642294431288815669E-3, b0, b1, b2);
        besselm1nextcheb(y, -1.05640848946261981558E-2, b0, b1, b2);
        besselm1nextcheb(y, 2.47264490306265168283E-2, b0, b1, b2);
        besselm1nextcheb(y, -5.29459812080949914269E-2, b0, b1, b2);
        besselm1nextcheb(y, 1.02643658689847095384E-1, b0, b1, b2);
        besselm1nextcheb(y, -1.76416518357834055153E-1, b0, b1, b2);
        besselm1nextcheb(y, 2.52587186443633654823E-1, b0, b1, b2);
        v = 0.5*(b0-b2);
        z = v*z*exp(z);
    }
    else
    {
        y = 32.0/z-2.0;
        besselm1firstcheb(7.51729631084210481353E-18, b0, b1, b2);
        besselm1nextcheb(y, 4.41434832307170791151E-18, b0, b1, b2);
        besselm1nextcheb(y, -4.65030536848935832153E-17, b0, b1, b2);
        besselm1nextcheb(y, -3.20952592199342395980E-17, b0, b1, b2);
        besselm1nextcheb(y, 2.96262899764595013876E-16, b0, b1, b2);
        besselm1nextcheb(y, 3.30820231092092828324E-16, b0, b1, b2);
        besselm1nextcheb(y, -1.88035477551078244854E-15, b0, b1, b2);
        besselm1nextcheb(y, -3.81440307243700780478E-15, b0, b1, b2);
        besselm1nextcheb(y, 1.04202769841288027642E-14, b0, b1, b2);
        besselm1nextcheb(y, 4.27244001671195135429E-14, b0, b1, b2);
        besselm1nextcheb(y, -2.10154184277266431302E-14, b0, b1, b2);
        besselm1nextcheb(y, -4.08355111109219731823E-13, b0, b1, b2);
        besselm1nextcheb(y, -7.19855177624590851209E-13, b0, b1, b2);
        besselm1nextcheb(y, 2.03562854414708950722E-12, b0, b1, b2);
        besselm1nextcheb(y, 1.41258074366137813316E-11, b0, b1, b2);
        besselm1nextcheb(y, 3.25260358301548823856E-11, b0, b1, b2);
        besselm1nextcheb(y, -1.89749581235054123450E-11, b0, b1, b2);
        besselm1nextcheb(y, -5.58974346219658380687E-10, b0, b1, b2);
        besselm1nextcheb(y, -3.83538038596423702205E-9, b0, b1, b2);
        besselm1nextcheb(y, -2.63146884688951950684E-8, b0, b1, b2);
        besselm1nextcheb(y, -2.51223623787020892529E-7, b0, b1, b2);
        besselm1nextcheb(y, -3.88256480887769039346E-6, b0, b1, b2);
        besselm1nextcheb(y, -1.10588938762623716291E-4, b0, b1, b2);
        besselm1nextcheb(y, -9.76109749136146840777E-3, b0, b1, b2);
        besselm1nextcheb(y, 7.78576235018280120474E-1, b0, b1, b2);
        v = 0.5*(b0-b2);
        z = v*exp(z)/sqrt(z);
    }
    if( x<0 )
    {
        z = -z;
    }
    result = z;
    return result;
}


/*************************************************************************
Modified Bessel function, second kind, order zero

Returns modified Bessel function of the second kind
of order zero of the argument.

The range is partitioned into the two intervals [0,8] and
(8, infinity).  Chebyshev polynomial expansions are employed
in each interval.

ACCURACY:

Tested at 2000 random points between 0 and 8.  Peak absolute
error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      0, 30       30000       1.2e-15     1.6e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
double besselk0(double x)
{
    double result;
    double y;
    double z;
    double v;
    double b0;
    double b1;
    double b2;

    ap::ap_error::make_assertion(x>0, "Domain error in BesselK0: x<=0");
    if( x<=2 )
    {
        y = x*x-2.0;
        besselmfirstcheb(1.37446543561352307156E-16, b0, b1, b2);
        besselmnextcheb(y, 4.25981614279661018399E-14, b0, b1, b2);
        besselmnextcheb(y, 1.03496952576338420167E-11, b0, b1, b2);
        besselmnextcheb(y, 1.90451637722020886025E-9, b0, b1, b2);
        besselmnextcheb(y, 2.53479107902614945675E-7, b0, b1, b2);
        besselmnextcheb(y, 2.28621210311945178607E-5, b0, b1, b2);
        besselmnextcheb(y, 1.26461541144692592338E-3, b0, b1, b2);
        besselmnextcheb(y, 3.59799365153615016266E-2, b0, b1, b2);
        besselmnextcheb(y, 3.44289899924628486886E-1, b0, b1, b2);
        besselmnextcheb(y, -5.35327393233902768720E-1, b0, b1, b2);
        v = 0.5*(b0-b2);
        v = v-log(0.5*x)*besseli0(x);
    }
    else
    {
        z = 8.0/x-2.0;
        besselmfirstcheb(5.30043377268626276149E-18, b0, b1, b2);
        besselmnextcheb(z, -1.64758043015242134646E-17, b0, b1, b2);
        besselmnextcheb(z, 5.21039150503902756861E-17, b0, b1, b2);
        besselmnextcheb(z, -1.67823109680541210385E-16, b0, b1, b2);
        besselmnextcheb(z, 5.51205597852431940784E-16, b0, b1, b2);
        besselmnextcheb(z, -1.84859337734377901440E-15, b0, b1, b2);
        besselmnextcheb(z, 6.34007647740507060557E-15, b0, b1, b2);
        besselmnextcheb(z, -2.22751332699166985548E-14, b0, b1, b2);
        besselmnextcheb(z, 8.03289077536357521100E-14, b0, b1, b2);
        besselmnextcheb(z, -2.98009692317273043925E-13, b0, b1, b2);
        besselmnextcheb(z, 1.14034058820847496303E-12, b0, b1, b2);
        besselmnextcheb(z, -4.51459788337394416547E-12, b0, b1, b2);
        besselmnextcheb(z, 1.85594911495471785253E-11, b0, b1, b2);
        besselmnextcheb(z, -7.95748924447710747776E-11, b0, b1, b2);
        besselmnextcheb(z, 3.57739728140030116597E-10, b0, b1, b2);
        besselmnextcheb(z, -1.69753450938905987466E-9, b0, b1, b2);
        besselmnextcheb(z, 8.57403401741422608519E-9, b0, b1, b2);
        besselmnextcheb(z, -4.66048989768794782956E-8, b0, b1, b2);
        besselmnextcheb(z, 2.76681363944501510342E-7, b0, b1, b2);
        besselmnextcheb(z, -1.83175552271911948767E-6, b0, b1, b2);
        besselmnextcheb(z, 1.39498137188764993662E-5, b0, b1, b2);
        besselmnextcheb(z, -1.28495495816278026384E-4, b0, b1, b2);
        besselmnextcheb(z, 1.56988388573005337491E-3, b0, b1, b2);
        besselmnextcheb(z, -3.14481013119645005427E-2, b0, b1, b2);
        besselmnextcheb(z, 2.44030308206595545468E0, b0, b1, b2);
        v = 0.5*(b0-b2);
        v = v*exp(-x)/sqrt(x);
    }
    result = v;
    return result;
}


/*************************************************************************
Modified Bessel function, second kind, order one

Computes the modified Bessel function of the second kind
of order one of the argument.

The range is partitioned into the two intervals [0,2] and
(2, infinity).  Chebyshev polynomial expansions are employed
in each interval.

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      0, 30       30000       1.2e-15     1.6e-16

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
double besselk1(double x)
{
    double result;
    double y;
    double z;
    double v;
    double b0;
    double b1;
    double b2;

    z = 0.5*x;
    ap::ap_error::make_assertion(z>0, "Domain error in K1");
    if( x<=2 )
    {
        y = x*x-2.0;
        besselm1firstcheb(-7.02386347938628759343E-18, b0, b1, b2);
        besselm1nextcheb(y, -2.42744985051936593393E-15, b0, b1, b2);
        besselm1nextcheb(y, -6.66690169419932900609E-13, b0, b1, b2);
        besselm1nextcheb(y, -1.41148839263352776110E-10, b0, b1, b2);
        besselm1nextcheb(y, -2.21338763073472585583E-8, b0, b1, b2);
        besselm1nextcheb(y, -2.43340614156596823496E-6, b0, b1, b2);
        besselm1nextcheb(y, -1.73028895751305206302E-4, b0, b1, b2);
        besselm1nextcheb(y, -6.97572385963986435018E-3, b0, b1, b2);
        besselm1nextcheb(y, -1.22611180822657148235E-1, b0, b1, b2);
        besselm1nextcheb(y, -3.53155960776544875667E-1, b0, b1, b2);
        besselm1nextcheb(y, 1.52530022733894777053E0, b0, b1, b2);
        v = 0.5*(b0-b2);
        result = log(z)*besseli1(x)+v/x;
    }
    else
    {
        y = 8.0/x-2.0;
        besselm1firstcheb(-5.75674448366501715755E-18, b0, b1, b2);
        besselm1nextcheb(y, 1.79405087314755922667E-17, b0, b1, b2);
        besselm1nextcheb(y, -5.68946255844285935196E-17, b0, b1, b2);
        besselm1nextcheb(y, 1.83809354436663880070E-16, b0, b1, b2);
        besselm1nextcheb(y, -6.05704724837331885336E-16, b0, b1, b2);
        besselm1nextcheb(y, 2.03870316562433424052E-15, b0, b1, b2);
        besselm1nextcheb(y, -7.01983709041831346144E-15, b0, b1, b2);
        besselm1nextcheb(y, 2.47715442448130437068E-14, b0, b1, b2);
        besselm1nextcheb(y, -8.97670518232499435011E-14, b0, b1, b2);
        besselm1nextcheb(y, 3.34841966607842919884E-13, b0, b1, b2);
        besselm1nextcheb(y, -1.28917396095102890680E-12, b0, b1, b2);
        besselm1nextcheb(y, 5.13963967348173025100E-12, b0, b1, b2);
        besselm1nextcheb(y, -2.12996783842756842877E-11, b0, b1, b2);
        besselm1nextcheb(y, 9.21831518760500529508E-11, b0, b1, b2);
        besselm1nextcheb(y, -4.19035475934189648750E-10, b0, b1, b2);
        besselm1nextcheb(y, 2.01504975519703286596E-9, b0, b1, b2);
        besselm1nextcheb(y, -1.03457624656780970260E-8, b0, b1, b2);
        besselm1nextcheb(y, 5.74108412545004946722E-8, b0, b1, b2);
        besselm1nextcheb(y, -3.50196060308781257119E-7, b0, b1, b2);
        besselm1nextcheb(y, 2.40648494783721712015E-6, b0, b1, b2);
        besselm1nextcheb(y, -1.93619797416608296024E-5, b0, b1, b2);
        besselm1nextcheb(y, 1.95215518471351631108E-4, b0, b1, b2);
        besselm1nextcheb(y, -2.85781685962277938680E-3, b0, b1, b2);
        besselm1nextcheb(y, 1.03923736576817238437E-1, b0, b1, b2);
        besselm1nextcheb(y, 2.72062619048444266945E0, b0, b1, b2);
        v = 0.5*(b0-b2);
        result = exp(-x)*v/sqrt(x);
    }
    return result;
}


/*************************************************************************
Modified Bessel function, second kind, integer order

Returns modified Bessel function of the second kind
of order n of the argument.

The range is partitioned into the two intervals [0,9.55] and
(9.55, infinity).  An ascending power series is used in the
low range, and an asymptotic expansion in the high range.

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      0,30        90000       1.8e-8      3.0e-10

Error is high only near the crossover point x = 9.55
between the two expansions used.

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
*************************************************************************/
double besselkn(int nn, double x)
{
    double result;
    double k;
    double kf;
    double nk1f;
    double nkf;
    double zn;
    double t;
    double s;
    double z0;
    double z;
    double ans;
    double fn;
    double pn;
    double pk;
    double zmn;
    double tlg;
    double tox;
    int i;
    int n;
    double eul;

    eul = 5.772156649015328606065e-1;
    if( nn<0 )
    {
        n = -nn;
    }
    else
    {
        n = nn;
    }
    ap::ap_error::make_assertion(n<=31, "Overflow in BesselKN");
    ap::ap_error::make_assertion(x>0, "Domain error in BesselKN");
    if( x<=9.55 )
    {
        ans = 0.0;
        z0 = 0.25*x*x;
        fn = 1.0;
        pn = 0.0;
        zmn = 1.0;
        tox = 2.0/x;
        if( n>0 )
        {
            pn = -eul;
            k = 1.0;
            for(i = 1; i <= n-1; i++)
            {
                pn = pn+1.0/k;
                k = k+1.0;
                fn = fn*k;
            }
            zmn = tox;
            if( n==1 )
            {
                ans = 1.0/x;
            }
            else
            {
                nk1f = fn/n;
                kf = 1.0;
                s = nk1f;
                z = -z0;
                zn = 1.0;
                for(i = 1; i <= n-1; i++)
                {
                    nk1f = nk1f/(n-i);
                    kf = kf*i;
                    zn = zn*z;
                    t = nk1f*zn/kf;
                    s = s+t;
                    ap::ap_error::make_assertion(ap::maxrealnumber-fabs(t)>fabs(s), "Overflow in BesselKN");
                    ap::ap_error::make_assertion(!(tox>1.0&&ap::maxrealnumber/tox<zmn), "Overflow in BesselKN");
                    zmn = zmn*tox;
                }
                s = s*0.5;
                t = fabs(s);
                ap::ap_error::make_assertion(!(zmn>1.0&&ap::maxrealnumber/zmn<t), "Overflow in BesselKN");
                ap::ap_error::make_assertion(!(t>1.0&&ap::maxrealnumber/t<zmn), "Overflow in BesselKN");
                ans = s*zmn;
            }
        }
        tlg = 2.0*log(0.5*x);
        pk = -eul;
        if( n==0 )
        {
            pn = pk;
            t = 1.0;
        }
        else
        {
            pn = pn+1.0/n;
            t = 1.0/fn;
        }
        s = (pk+pn-tlg)*t;
        k = 1.0;
        do
        {
            t = t*(z0/(k*(k+n)));
            pk = pk+1.0/k;
            pn = pn+1.0/(k+n);
            s = s+(pk+pn-tlg)*t;
            k = k+1.0;
        }
        while(fabs(t/s)>ap::machineepsilon);
        s = 0.5*s/zmn;
        if( n%2!=0 )
        {
            s = -s;
        }
        ans = ans+s;
        result = ans;
        return result;
    }
    if( x>log(ap::maxrealnumber) )
    {
        result = 0;
        return result;
    }
    k = n;
    pn = 4.0*k*k;
    pk = 1.0;
    z0 = 8.0*x;
    fn = 1.0;
    t = 1.0;
    s = t;
    nkf = ap::maxrealnumber;
    i = 0;
    do
    {
        z = pn-pk*pk;
        t = t*z/(fn*z0);
        nk1f = fabs(t);
        if( i>=n&&nk1f>nkf )
        {
            break;
        }
        nkf = nk1f;
        s = s+t;
        fn = fn+1.0;
        pk = pk+2.0;
        i = i+1;
    }
    while(fabs(t/s)>ap::machineepsilon);
    result = exp(-x)*sqrt(ap::pi()/(2.0*x))*s;
    return result;
}


/*************************************************************************
Internal subroutine

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
static void besselmfirstcheb(double c, double& b0, double& b1, double& b2)
{

    b0 = c;
    b1 = 0.0;
    b2 = 0.0;
}


/*************************************************************************
Internal subroutine

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
static void besselmnextcheb(double x,
     double c,
     double& b0,
     double& b1,
     double& b2)
{

    b2 = b1;
    b1 = b0;
    b0 = x*b1-b2+c;
}


/*************************************************************************
Internal subroutine

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
static void besselm1firstcheb(double c, double& b0, double& b1, double& b2)
{

    b0 = c;
    b1 = 0.0;
    b2 = 0.0;
}


/*************************************************************************
Internal subroutine

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
static void besselm1nextcheb(double x,
     double c,
     double& b0,
     double& b1,
     double& b2)
{

    b2 = b1;
    b1 = b0;
    b0 = x*b1-b2+c;
}


static void besselasympt0(double x, double& pzero, double& qzero)
{
    double xsq;
    double p2;
    double q2;
    double p3;
    double q3;

    xsq = 64.0/(x*x);
    p2 = 0.0;
    p2 = 2485.271928957404011288128951+xsq*p2;
    p2 = 153982.6532623911470917825993+xsq*p2;
    p2 = 2016135.283049983642487182349+xsq*p2;
    p2 = 8413041.456550439208464315611+xsq*p2;
    p2 = 12332384.76817638145232406055+xsq*p2;
    p2 = 5393485.083869438325262122897+xsq*p2;
    q2 = 1.0;
    q2 = 2615.700736920839685159081813+xsq*q2;
    q2 = 156001.7276940030940592769933+xsq*q2;
    q2 = 2025066.801570134013891035236+xsq*q2;
    q2 = 8426449.050629797331554404810+xsq*q2;
    q2 = 12338310.22786324960844856182+xsq*q2;
    q2 = 5393485.083869438325560444960+xsq*q2;
    p3 = -0.0;
    p3 = -4.887199395841261531199129300+xsq*p3;
    p3 = -226.2630641933704113967255053+xsq*p3;
    p3 = -2365.956170779108192723612816+xsq*p3;
    p3 = -8239.066313485606568803548860+xsq*p3;
    p3 = -10381.41698748464093880530341+xsq*p3;
    p3 = -3984.617357595222463506790588+xsq*p3;
    q3 = 1.0;
    q3 = 408.7714673983499223402830260+xsq*q3;
    q3 = 15704.89191515395519392882766+xsq*q3;
    q3 = 156021.3206679291652539287109+xsq*q3;
    q3 = 533291.3634216897168722255057+xsq*q3;
    q3 = 666745.4239319826986004038103+xsq*q3;
    q3 = 255015.5108860942382983170882+xsq*q3;
    pzero = p2/q2;
    qzero = 8*p3/q3/x;
}


static void besselasympt1(double x, double& pzero, double& qzero)
{
    double xsq;
    double p2;
    double q2;
    double p3;
    double q3;

    xsq = 64.0/(x*x);
    p2 = -1611.616644324610116477412898;
    p2 = -109824.0554345934672737413139+xsq*p2;
    p2 = -1523529.351181137383255105722+xsq*p2;
    p2 = -6603373.248364939109255245434+xsq*p2;
    p2 = -9942246.505077641195658377899+xsq*p2;
    p2 = -4435757.816794127857114720794+xsq*p2;
    q2 = 1.0;
    q2 = -1455.009440190496182453565068+xsq*q2;
    q2 = -107263.8599110382011903063867+xsq*q2;
    q2 = -1511809.506634160881644546358+xsq*q2;
    q2 = -6585339.479723087072826915069+xsq*q2;
    q2 = -9934124.389934585658967556309+xsq*q2;
    q2 = -4435757.816794127856828016962+xsq*q2;
    p3 = 35.26513384663603218592175580;
    p3 = 1706.375429020768002061283546+xsq*p3;
    p3 = 18494.26287322386679652009819+xsq*p3;
    p3 = 66178.83658127083517939992166+xsq*p3;
    p3 = 85145.16067533570196555001171+xsq*p3;
    p3 = 33220.91340985722351859704442+xsq*p3;
    q3 = 1.0;
    q3 = 863.8367769604990967475517183+xsq*q3;
    q3 = 37890.22974577220264142952256+xsq*q3;
    q3 = 400294.4358226697511708610813+xsq*q3;
    q3 = 1419460.669603720892855755253+xsq*q3;
    q3 = 1819458.042243997298924553839+xsq*q3;
    q3 = 708712.8194102874357377502472+xsq*q3;
    pzero = p2/q2;
    qzero = 8*p3/q3/x;
}



